clear
%% Problem parameters
%   'EqualMinima','Himmelblau','Sixhump','ModifiedRastrigin','Vincent',
%   'Shubert','Composition','IncreasingMinima','Rastrigin','Schaffer'.
ProblemSet.FuncType = 'IncreasingMinima';
ProblemSet.dimension = 1;
% ProblemSet.k = [3,3,3];
% ProblemSet.BaseFunc = 3;

[d, x_bound, ObjFunc, OptSet, precision_init, precision, ~, MaximalBudget, FileName]...
    = problem_setting(ProblemSet);

MaximalBudget = 1000;
iepsilon = 4;
epsilon = [0.1,0.01,0.001,0.0001];
OptSet.epsilon = epsilon(iepsilon);
%% problem
NewRegionNum = 2;
MaxPartitionDepth_d = ceil(log((x_bound(2,:)-x_bound(1,:))./precision_init)...
    ./log(NewRegionNum));
MaxPartitionDepth = sum(MaxPartitionDepth_d)
ClearRadius0 = (x_bound(2,:)-x_bound(1,:))./(NewRegionNum.^(MaxPartitionDepth_d));
ClearRadius = 2*min(ClearRadius0);

Problem.Dimension = d;
Problem.Domain = reshape(x_bound,[1,d*2]);
Problem.Partition = @(Region,SampleSet)Partition_BoxConstraint(Region, NewRegionNum, MaxPartitionDepth, SampleSet);
Problem.Sampling = @(n, region)Sampling_BoxConstraint(n, region, ObjFunc);

%% optimization
% Algorithm parameters
AlgorithmM.n0 = 4;
AlgorithmM.QuantileLevel = 0.3;
AlgorithmM.NewBudget = 3;
AlgorithmM.MaxSampleSize = 10; %5*(d+1);
AlgorithmM.radius = ClearRadius;
AlgorithmM.local_search = @(starts,step)LocalSearch_BoxConstraint(ObjFunc,x_bound,starts,step,precision(iepsilon)/2,MaximalBudget);
AlgorithmM.StopCriteria = [1,MaximalBudget];
% AlgorithmM.StopCriteria = [2,size(OptSet.sol,1)];
% AlgorithmM.StopCriteria = [3,5000];
%% OPTIMIZATION
rep = 500;
OptObtained = cell(rep,1);
TotalBudget1 = zeros(rep,1);
TotalBudget2 = zeros(rep,1);
BudgetOptima = zeros(rep,size(OptSet.sol,1));
Times = zeros(rep,1);
%%
for i=1:rep
    disp(['replication: ',num2str(i),'.'])
    tic;
    [optima, ~, SampleSet, ~, ls_budget, BudgetForEachOptima] = prsmmo_ls( Problem, AlgorithmM, OptSet );
    TotalBudget1(i) = size(SampleSet,1) - ls_budget;
    TotalBudget2(i) = ls_budget;
    Times(i) = toc;
    OptObtained{i} = optima;
    BudgetOptima(i,:) = BudgetForEachOptima';
    disp(['time used: ',num2str(Times(i)),'.'])
end

%% save data
clear optima SampleSet ls_budget
save(FileName)

%% FIGURES
NumOptPlot = size(OptSet.sol,1);  % size(OptSet.sol,1);
BudgetOptima(BudgetOptima==0)=inf;
% opt_i_budget = zeros(rep,size(NumOptPlot,1));  % the budget required for each optimum
% for i=1:rep
%     for j = 1:NumOptPlot
%         dis = sqrt(sum((OptObtained{i}(:,3:end)-repmat(OptSet.sol(j,:),[size(OptObtained{i},1),1])).^2,2));
%         id = find(dis<OptSet.radius);
%         if isempty(id)
%             opt_i_budget(i,j) = Inf;
%         else
%             [~,id_min] = min(OptObtained{i}(id,1));
%             id = id(id_min);
%             opt_i_budget(i,j) = OptObtained{i}(id,1);
%         end
%         if ~isempty(id) && OptObtained{i}(id,2)>(OptSet.fval(j)+epsilon(iepsilon))
%             error('the optimum is not found.')
%         end
%     end
% end
%% seperate all optima
% figure()
% hold on
% legend_i = cell(1,NumOptPlot);
% for i = 1:NumOptPlot
%     N = 1:MaximalBudget;
%     NumOpt = zeros(1,length(N));
%     for j=1:length(N)
%         NumOpt(j) = sum(BudgetOptima(:,i)<N(j));
%     end
%     plot(N,NumOpt./rep)
%     legend_i{i} = ['f(x',num2str(i),')=',num2str(OptSet.fval(i),'%2.2f')];
% end
% legend(legend_i)
% xlabel('Total budget')
% ylabel('Probability of being found')
% set(gca,'Fontname','Times New Roman','FontSize',16);
% hold off
%% combine optima with the same quality
OptSet.fval = round(OptSet.fval./epsilon(iepsilon)./10).*epsilon(iepsilon).*10;
DiffOpt = unique(OptSet.fval(1:NumOptPlot));
NumDiffOpt = size(DiffOpt,1);
figure()
hold on
legend_i = cell(1,NumDiffOpt);
for i = 1:NumDiffOpt
    N = 1:MaximalBudget;
    SameOpt = find(OptSet.fval(1:NumOptPlot)==DiffOpt(i));
    NumOpt = zeros(1,length(N));
    for j=1:length(N)
        NumOpt(j) = sum(sum(BudgetOptima(:,SameOpt)<N(j)));
    end
    NumOpt = NumOpt./length(SameOpt);
    plot(N,NumOpt./rep)
    legend_i{i} = ['f(x) = ',num2str(OptSet.fval(SameOpt(1)),'%2.2f'),' : ', num2str(length(SameOpt))];
end
legend(legend_i)
xlabel('Total budget size')
ylabel('Probability of being found')
set(gca,'Fontname','Times New Roman','FontSize',16);
hold off
